Accelerating the Solution of Linear Systems by Iterative Refinement in Three Precisions∗

We propose a general algorithm for solving a n×n nonsingular linear system Ax = b based on iterative refinement with three precisions. The working precision is combined with possibly different precisions for solving for the correction term and for computing the residuals. Via rounding error analysis of the algorithm we derive sufficient conditions for convergence and bounds for the attainable normwise forward error and normwise and componentwise backward errors. Our results generalize and unify many existing rounding error analyses for iterative refinement. With single precision as the working precision, we show that by using LU factorization in IEEE half precision as the solver and calculating the residuals in double precision it is possible to solve Ax = b to full single precision accuracy for condition numbers κ2(A) ≤ 104, with the O(n3) part of the computations carried out entirely in half precision. We show further that by solving the correction equations by GMRES preconditioned by the LU factors the restriction on the condition number can be weakened to κ2(A) ≤ 108, although in general there is no guarantee that GMRES will converge quickly. Taking for comparison a standard Ax = b solver that uses LU factorization in single precision, these results suggest that on architectures for which half precision is efficiently implemented it will be possible to solve certain linear systems Ax = b up to twice as fast and to greater accuracy. Analogous results are given with double precision as the working precision.